home *** CD-ROM | disk | FTP | other *** search
/ Nautilus 1992 July / Nautilus-3-8 / Nautilus-3-8.bin / Tools & Utilities / Techy Stuff / Source ƒ / sky source ƒ / NUTATE.C < prev    next >
Encoding:
C/C++ Source or Header  |  1992-06-16  |  2.5 KB  |  130 lines
  1. #include "sky.h"
  2.  
  3. extern struct nuttab
  4. {
  5.     float f[2];
  6.     char c[5];
  7. } nuttab[];
  8.  
  9. nutate()
  10. {
  11.  
  12. /*
  13.  *    uses radian, radsec, eday
  14.  *    sets psi, eps, dpsi, deps, obliq, gst, tobliq
  15.  */
  16.  
  17.     double msun, mnom, noded, dmoon;
  18.     struct nuttab *pp;
  19.     double arg;
  20.  
  21. /*
  22.  *    nutation of the equinoxes is a wobble of the pole
  23.  *    of the earths rotation whose magnitude is about
  24.  *    9 seconds of arc and whose period is about 18.6 years.
  25.  *
  26.  *    it depends upon the pull of the sun and moon on the
  27.  *    equatorial bulge of the earth.
  28.  *
  29.  *    psi and eps are the two angles which specify the
  30.  *    true pole with respect to the mean pole.
  31.  *
  32.  *    all coeffieients are from Exp. Supp. pp.44-45
  33.  */
  34.  
  35.     pp = &nuttab[0];
  36.     mnom = 296.104608 + 13.0649924465*eday + 9.192e-3*capt2
  37.          + 14.38e-6*capt3;
  38.     msun = 358.475833 + .9856002669*eday - .150e-3*capt2
  39.          - 3.33e-6*capt3;
  40.     noded = 11.250889 + 13.2293504490*eday - 3.211e-3*capt2
  41.          - 0.33e-6*capt3;
  42.     dmoon = 350.737486 + 12.1907491914*eday - 1.436e-3*capt2
  43.          + 1.89e-6*capt3;
  44.     node = 259.183275 - .0529539222*eday + 2.078e-3*capt2
  45.          + 2.22e-6*capt3;
  46.     mnom *= radian;
  47.     msun *= radian;
  48.     noded *= radian;
  49.     dmoon *= radian;
  50.     node *= radian;
  51.  
  52. /*
  53.  *    long period terms
  54.  */
  55.  
  56.  
  57.     psi = -(17.2327+.01737*capt)*sin(node);
  58.     eps = (9.2100+0.00091*capt)*cos(node);
  59.     for(;;){
  60.         if(pp->f[0]==0.){
  61.             pp++;
  62.             break;
  63.         }
  64.         arg = pp->c[0]*mnom + pp->c[1]*msun + pp->c[2]*noded
  65.             + pp->c[3]*dmoon + pp->c[4]*node;
  66.         psi += pp->f[0]*sin(arg);
  67.         eps += pp->f[1]*cos(arg);
  68.         pp++;
  69.     }
  70.  
  71. /*
  72.  *    short period terms
  73.  */
  74.  
  75.     dpsi = 0.;
  76.     deps = 0.;
  77.     if((flags&APPARENT)==0){
  78.         for(;;){
  79.             if(pp->f[0]==0.){
  80.                 pp++;
  81.                 break;
  82.             }
  83.             arg = pp->c[0]*mnom + pp->c[1]*msun + pp->c[2]*noded
  84.                 + pp->c[3]*dmoon + pp->c[4]*node;
  85.             dpsi += pp->f[0]*sin(arg);
  86.             deps += pp->f[1]*cos(arg);
  87.             pp++;
  88.         }
  89.     }
  90.  
  91.  
  92.     psi = (psi+dpsi)*radsec;
  93.     eps = (eps+deps)*radsec;
  94.     dpsi *= radsec;
  95.     deps *= radsec;
  96.  
  97.     obliq = 23.452294 - .0130125*capt - 1.64e-6*capt2
  98.          + 0.503e-6*capt3;
  99.     obliq *= radian;
  100.     tobliq = obliq + eps;
  101.  
  102.     gst = 99.690983 + 360.9856473354*eday + .000387*capt2;
  103.     gst -= 180.;
  104.     gst = fmod(gst, 360.);
  105.     if(gst < 0.)
  106.         gst += 360.;
  107.     gst *= radian;
  108.     gst += psi*cos(obliq);
  109.  
  110. /*
  111.  * print true obliquity, mean obliquity, mean GST, and true GST.
  112.  */
  113.  
  114. /*
  115.     printf("Mean. G. Sid. Time: ");
  116.     prhms(gst-psi*cos(obliq),3);
  117.     printf("\n");
  118.     printf("App. G. Sid. Time: ");
  119.     prhms(gst,3);
  120.     printf("\n");
  121.  
  122.     printf("Goobie: %6.3f %6.3f %6.3f %6.3f\n", psi/radsec, eps/radsec, dpsi/radsec, deps/radsec);
  123.  
  124.     printf("Obliquity: ");
  125.     prdms(tobliq,3);
  126.     printf("\n");
  127. */
  128.  
  129. }
  130.